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Abstract 

Using the optimized effective potential method in conjunction with 
the semi-analytical approximation due to Krieger, Li and lafrate, we 
have performed fully self-consistent exact exchange-only 
density-functional calculations for diatomic molecules with a fully 
numerical basis-set-free molecular code. The results are very similar 
to the ones obtained with the Hartree Fock approach. Furthermore 
we present results for ground states of positive atomic ions including 
correlation contributions in the approximation of Colle and Salvetti. It 
is found that the scheme performs significantly better than 
conventional Kohn-Sham calculations. 



1 Introduction 

Since its development by Talman and Shadwick [|l[], following the original idea of 
Sharp and Norton |^], the optimized effective potential (OEP) method has been 
recognized ||3|, ^ as the exact implementation of exchange-only density functional 
theory (DFT) [|5|, ^, ^]. Due to the rather large computational effort involved, 
this scheme has not been used extensively. In the last years, however, the situation 
has changed. Using an accurate analytical approximation due to Krieger, Li and 
lafrate (KLI) 1^, |l^, 11] the effort involved in numerical calculations based on the 



OEP has become comparable to conventional Kohn-Sham calculations, while the 
gain in accuracy is considerable. 

In the following, we will briefly review the theoretical foundations of the OEP 
method and the KLI approximation and present, in section ^, applications of the 
method to molecular systems. Finally, in section ^ calculations for atomic systems 



including correlation effects are discussed and compared with conventional Kohn- 
Sham results. 



2 Basic formalism 

We start from ordinary spin DFT [|I2|, |T3|], where the basic variables are the spin- 
up and spin-down densities /0|(r) and Pj(r), respectively. They are obtained by 
self-consistently solving the single-particle Schrodinger equations (atomic units are 
used throughout) 

( — ^ + ^a[/3T>/'i](r) 1 f/JiaW =eja'^if7(r) j = l,...,N^ =T, i (1) 

where 

p^{r)=Y^\^Ur)\', (2) 

i=l 

The Kohn-Sham potentials ^^(r) may be written in the usual way as 

K(r) = v,^,{r) + J d\' + Kea(r), (3) 

p(r) = P'^i^) (4) 

where fext(r) represents the Coulomb potential of the nuclei and T^co-lr) is a local 
exchange-correlation (xc) potential formally defined as functional derivative of the 
xc energy 

V.c.[r).- ^^^^^^ . (5) 

In order to understand the nature of the OEP method we recall that the Hohenberg- 
Kohn theorem |P, applied to non-interacting systems, ensures that the ground- 
state determinant and hence all occupied orbitals are unique functionals of the spin 
densities: 

iPja{r) = ipja[p^,Pl]{r). (6) 

As a consequence of (^), any orbital functional Exdi^jr}] is an implicit func- 
tional of p| and p|, provided that the orbitals come from a local potential. 
The starting point of the so-called OEP method is the total energy functional 

+ / d^r t;ext(r)p(r) 



+ '-[d\[ d\' 



|r — r' 

nOEP 



[Wn)] (7) 
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where, in contrast to ordinary spin DFT, the xc energy is an explicit (approximate) 
functional of spin orbitals and therefore only an implicit functional of the spin 
densities and pj. As a consequence of this fact, the calculation of the xc 
potentials from Eq. (^) is somewhat more complicated: We use the chain rule for 
functional derivatives to obtain 

and, by applying the functional chain rule once more, 

(9) 

The last term on the right-hand side is readily identified with the inverse Xs ^ (''j 
of the density response function of a system of non-interacting particles 

This quantity is diagonal with respect to the spin variables so that Eq. (^) reduces 
to 



(11) 

Acting with the response operator ([lO| ) on both sides of Eq. (|Il| ) one obtains 

(12) 

Finally, the second functional derivative on the right-hand side of Eq. (|l^) is 
calculated using first-order perturbation theory. This yields 

Using the fact that the Kohn-Sham response function can be written as 

Xs. (r, r') = E E + C.C. (14) 



the integral equation (12) takes the form 



i=l ■ 



E / d'r h/°r(r') - w(r') G,,,, (r', r) ^^Ar)^l{r') + c.c. = (15) 
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where 



and 



oo 



(17) 

fc=i ^«<^ ^'^'^ 

The derivation of the OEP integral equation (p^) described here was first given 



by Gorling and Levy It is important to note that the same expression results 
[||, ^, ^ if one demands that the local one-particle potential appearing in Eq. (||) 
be the optimized one yielding orbitals minimizing the total energy functional (0), 
i.e. that 

5E?r 



= 0. (18) 

y^yOEP 

The main advantage of the OEP method is that it allows for the exact treatment 
of the exchange energy. Splitting up the total xc-functional into an exchange and 
a correlation part 

Eg^'' iWjr}] = + E, (19) 

we can use the exact Fock expression 

tT=t,J. i,k=l ' ' 

Performing the functional derivative with respect to the orbitals one obtaines for 
the x-part nxio-(r) of the function u^ciai'^)'- 



f A^' ^&ippfl. (21) 



1 

The use of the exact exchange energy has several advantages over the con- 
ventional explicitly density dependent xc functional. Most importantly it ensures 
the correct — 1/r decay of the xc-potential for large r, reflecting the fact that it is 
self-interaction free. One has to emphasize that the OEP has the correct —1/r tail 
for all orbitals, i.e. for both the occupied and the unoccupied ones. By contrast, 
the conventional Hartree-Fock (HF) approach, which uses the same expression (^) 
for the exchange energy but a nonlocal potential defined via the equation 

(K<7 )ir) = -2^J d r , (22) 

j=i I I 

is self-interaction free only for the occupied orbitals. However, x-only OEP calcula- 
tions lH, |l^, 17, IS, |l9|, ^, i.e. with the approximation Ec [{^jr}] = 0, performed 



on atomic systems have shown that the results for various physical quantities of 
interest such as total ground-state energies are very similar to HF results despite 
the different nature of the corresponding x-potentials. As - by construction - the 
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HF scheme gives the variationally best, i.e. lowest, total energy, the x-only OEP 
solutions are always slightly higher in energy. 

The solution of the full integral equation (15) is numerically very demanding 



and has been achieved so far only for systems with spherical symmetry |16|, |17 



1£, IS, S]. Therefore, one has to resort to further approximations for practical 



reasons. Krieger, Li and lafrate have suggested the analytical approximation 

Gs,,. (r, r') « ^ (Mr - r') - ^i.(r)^:^(r')) (23) 



for the Green's-function-type quantity (|r7|). Substituting this into the integral 
equation ([T5|) and performing some algebra one arrives at the approximate equation 



V^^\r) = -i— Y: \^^Ar)\' [w(r) + (v^j^ - w) + c.c.l (24) 

where Uxcja denotes the average value of nxcjcr(r) taken over the density of the ja 
orbital, i.e. 

u^cja = / d^r \ipj^{r)\'^u^cja{r) (25) 



and similarly for V^^:^. In contrast to the exact integral equation ( p^ ) the KLI 
equation (|24l) can be solved explicitly for Vxca by multiplication with Pia{r) and 
subsequent integration. This leads to linear {N„ — 1) x {Na — 1) equations for the 
unknown constants {Vxda — 



with j = 1, . . . , A'^^ - 1 



(26) 



' Pair) 



and 



^^.^(r) := / d'r ^^^^4^ l¥'-(r)p| (w(r) + <,,.(r)) . (28) 
J Pair) fr[ 2 

The orbitals corresponding to the highest single-particle energy eigenvalues e^a 
have to be excluded from the linear equations (|26|) in order to ensure the correct 
long-range behaviour of 14°^^ (r) [^]. It is an important property of the KLI approx- 
imation that it is exact for two-particle systems, where one has only one electron 
per spin projection. In this case, the OEP integral equation (15) may be solved 
exactly to yield (^4|). Furthermore, for these systems the OEP is also identical with 



the HF potential 

At first sight, the KLI approximation ( ^3|) might appear rather crude. The 
final result ( |24l ) for the KLI potential, however, can also be understood as a 
well-defined mean-field approximation. Explicit calculations on atoms performed in 
the x-only limit [||, |l^, show that the KLI-approximation yields excellent results 
which differ only by a few ppm from the much more time-consuming exact solutions 
of the full integral equation (15). 
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3 Exchange-only results for molecular systems 



In order to demonstrate the validity of the KLI-approach for more complex systems, 
we have performed x-only OEP calculations for diatomic molecules in KLI approx- 
imation employing the exact exchange energy functional as defined by equation 
(^) and neglecting correlation effects. This approach will, in the following, be 
referred to as x-only KLI. Our calculations have been performed with a fully nu- 
merical basis-set-free code, developed from the Xa program written by Laaksonen, 



Sundholm and Pyykko |21, |2^, |2^. The code solves the one-particle Schrodinger 



equation for diatomic molecules 

" \R^\ ~ \B^\ + ^^(^^ + ^^f-'i^)) ^jA^) = ^^<r^U^l (29) 

where Ri denotes the location and Zi the nuclear charge of the z-th nucleus in the 
molecule. The partial differential equation is solved in prolate spheroidal coordinates 
on a two-dimensional mesh by a relaxation method, while the third variable, the 
azimuthal angle, is treated analytically. The Hartree potential 

VH(r) = / d\' (30) 
J |r — r'l 

and the functions nxicr(r) (cf. Eq. (21)) needed for the calculation of the exchange 



potential V^^{y) (cf. Eq. (|24|) ) are computed as solutions of a Poisson and Poisson- 
like equation, respectively. In this step, the same relaxation technique as for the 
solution of the Schrodinger equation ( ^9D is employed. Starting with an initial guess 
for the wave functions (/^^^(r), equations (|29|) , (|30|) , (|2l| ) together with (^) are 
iterated until self-consistent. A very detailed description of the code is given in 

In order to test the accuracy of the program, we have performed calculations 
on the Beryllium and Neon atom which are compared in Table || to exact results 
obtained with a one-dimensional atomic code. The results agree to all decimals 
given. Furthermore, from Table ^ it is evident that our program gives results for 
the two-electron molecules He2 and HeH+ which are identical to the HF ones 



obtained by Laaksonen et al ||2^ as they should be, as the HF and x-only KLI 
schemes are identical for these systems. 

For comparison, we have performed additional x-only calculations with two other 
approximations of T4cr(r) and E^, respectively. The first one of these, denoted by 
Slater in the following, uses - like the HF and x-only KLI method - the exact orbital 
representation of given in equation (^) but the averaged exchange potential 
due to Slater |25] given by 



= E ^U^)^M j d\' ^^jj^^, (31) 



which may be obtained from (|24|) by setting the constants Vxda — Uxcia equal to 
zero for all i. The other is the well known x-only local density approximation (LDA) 
of conventional DFT. As for the KLI calculations, we have successfully tested our 
implementations on atomic systems. 
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Beryllium 


Neon 




ID 


2D 


ID 


2D 


Grid 


400 


153 X 249 


400 


153 X 249 


Etot 


-14.5723 


-14.5723 


-128.5448 


-128.5448 


^HOMO 


-0.3089 


-0.3089 


-0.8494 


-0.8494 


< 1/r > 


2.1039 


2.1039 


3.1100 


3.1100 


< > 


4.3255 


4.3255 


0.9367 


0.9367 



Table 1: X-only KLI results for the Beryllium and Neon atom. ID denotes 
exact values obtained with our one-dimensional code, 2D the results from our 
two-dimensional code with one nuclear charge set equal to zero. All numbers 
in atomic units. 





H2 


HeH+ 


R 


1.4 


1.455 




HF 


x-only KLI 


HF 


x-only KLI 


Etot 


-1.133630 


-1.133630 


-2.933103 


-2.933103 




-0.594659 


-0.594659 


-1.637451 


-1.637451 


Q! 








-0.494460 


-0.494460 




0.243289 


0.243289 


0.373727 


0.373727 


QI 








-0.231525 


-0.231525 


QI 


0.090721 


0.090721 


0.173962 


0.173962 


< > 


2.573930 


2.573930 


1.340832 


1.340832 



Table 2: Comparison of results for two electron molecules with assumed bond 



length R. HF values from |2^. Qf, Q|, Q3, QI, denote the electronic con- 
tributions to the dipole, quadrupole, octopole and hexadecapole moments, re- 
spectively, determined from the molecular midpoint. All numbers in atomic 
units. 





HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-7.9874 


-7.9868 


-7.9811 


-7.7043 




-2.4452 


-2.0786 


-2.3977 


-1.7786 




-0.3017 


-0.3011 


-0.3150 


-0.1284 


Q! 


0.6531 


0.6440 


0.8614 


0.8679 


QI 


7.1282 


7.1365 


6.9657 


6.7717 


Qi 


2.9096 


2.9293 


3.0799 


2.6924 


QI 


16.0276 


16.1311 


15.5881 


15.0789 



Table 3: X-only results for LiH. HF values for bond length of 3.015 a.u. from 
1 24]. Present calculations performed on a 153 x 193 grid with bond length of 
3.015 a.u. All numbers in atomic units. 
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HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-25.1316 


-25.1290 


-25.1072 


-24.6299 




-7.6863 


-6.8624 


-7.4837 


-6.4715 




-0.6482 


-0.5856 


-0.6358 


-0.3956 


^Scr 


-0.3484 


-0.3462 


-0.3721 


-0.1626 


Q! 


5.3525 


5.3498 


5.2991 


5.3154 


Qi 


12.1862 


12.1416 


11.4720 


11.9542 


Qi 


15.6411 


15.5618 


14.3328 


14.0904 


QI 


25.8492 


25.4188 


25.2152 


21.9134 



Table 4: X-only results for BH. HF values for bond length of 2.336 a.u. from 
1 24]. Present calculations performed on a 193 x 265 grid with bond length of 
2.336 a.u. All numbers in atomic units. 





HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-100.0708 


-100.0675 


-100.0225 


-99.1512 


Elcr 


-26.2946 


-24.5116 


-25.6625 


-24.0209 


£2(7 


-1.6010 


-1.3994 


-1.4327 


-1.0448 


EScr 


-0.7682 


-0.7772 


-0.8167 


-0.4483 


EItt 


-0.6504 


-0.6453 


-0.6897 


-0.3109 


Qft 


-0.7561 


-0.8217 


-0.8502 


-0.6962 




1.7321 


1.8013 


1.8472 


1.7124 


Q!j°t 


-2.5924 


-2.7222 


-2.8782 


-2.4662 


Qf* 


5.0188 


5.1825 


5.3720 


4.7068 


< l/r >H 


6.1130 


6.0878 


6.0909 


6.0901 


< 1/r >F 


27.1682 


27.1622 


27.6049 


27.0289 



Table 5: X-only results for FH. HF values for bond length of 1.7328 a.u. from 
]24]. Present calculations performed on a 161 x 321 grid with bond length of 
1.7328 a.u. All numbers in atomic units. 





HF 


x-only KLI 


Slater 


X- 


only LDA 


Etot 


-5.72333 


-5.72332 


-5.72332 




-5.44740 


eio-g 


-0.92017 


-0.91929 


-0.91977 




-0.51970 




-0.91570 


-0.91566 


-0.91614 




-0.51452 


Qi 


31.36165 


31.35931 


31.35907 




31.35507 


QI 


245.8779 


245.8643 


245.8615 




245.8255 



Table 6: X-only results for He2. HF values for bond length of 5.6 a.u. from 
]24]. Present calculations performed on a 209 x 225 grid with bond length of 



5.6 a.u. All numbers in atomic units. 
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HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-14.8716 


-14.8706 


-14.8544 


-14.3970 


eio-g 


-2.4531 


-2.0276 


-2.3875 


-1.7869 


^Ictu 


-2.4528 


-2.0272 


-2.3873 


-1.7864 


£2<Tg 


-0.1820 


-0.1813 


-0.1989 


-0.0922 


Qi 


27.6362 


27.4993 


29.0014 


29.4401 


QI 


159.9924 


159.6809 


169.1300 


172.8505 



Table 7: X-only results for Li2. HF values for bond length of 5.051 a.u. from 
p^ . Present calculations performed on a 209 x 225 grid with bond length of 
5.051 a.u. All numbers in atomic units. 





HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-29.1337 


-29.1274 


-29.0939 


-28.4612 




-4.73150 


-4.09876 


-4.60353 


-3.78576 




-4.73147 


-4.09872 


-4.60351 


-3.87571 


^2<7g 


-0.39727 


-0.33452 


-0.37659 


-0.23067 




-0.24209 


-0.23489 


-0.26524 


-0.13163 


Qi 


46.0878 


46.2475 


43.4833 


46.2501 


QI 


261.774 


277.365 


249.135 


281.950 



Table 8: X-only results for Be2. HF values for bond length of 4.6 a.u. from 
1 24]. Present calculations performed on a 209 x 225 grid with bond length of 
4.6 a.u. All numbers in atomic units. 





HF 


x-only KLI 


Slater 


x-only LDA 


Etot 


-108.9936 


-108.9856 


-108.9110 


-107.7560 


Elcrg 


-15.6822 


-14.3722 


-15.2692 


-13.8950 


^Ictvl 


-15.6787 


-14.3709 


-15.2682 


-13.8936 


^2f7g 


-1.4726 


-1.3076 


-1.3316 


-0.9875 


e2o-u 


-0.7784 


-0.7453 


-0.7473 


-0.4434 


EScrg 


-0.6347 


-0.6305 


-0.6521 


-0.3335 


S^lvrg 


-0.6152 


-0.6818 


-0.6960 


-0.3887 




-0.9372 


-0.9488 


-1.1756 


-1.1643 




-7.3978 


-6.7476 


-7.1266 


-6.2553 


< l/r >N 


21.6543 


21.6439 


21.9749 


21.5820 



Table 9: X-only results for N2. HF values for bond length of 2.07 a.u. from 
1 24]. Present calculations performed on a 209 x 225 grid with bond length of 
2.07 a.u. All numbers in atomic units. 
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Results are given in Tables ^ through |^ for LiH, BH, FH, He2, \-\2, Be2 and 
N2. For each system we show the total ground state energy Etot. the various 
orbital energy eigenvalues e and the nonzero electronic contributions to the dipole, 
quadrupole, octopole and hexadecapole moments denoted by Qf , Ql, Q3 and Ql 
calculated from the geometrical center of the molecule, respectively, except for FH 
and N2, where the total moments (including nuclear contributions) calculated from 
the center of mass are given, denoted by Q^^*. For these two molecules we also 
present the expectation values of 1/r, denoted by < 1/r >, calculated at the nuclei. 

For all physical quantities of interest, i.e. for ExoT, the energies ehomo of the 
highest occupied orbitals and the multipole moments, the x-only KLI and HF results 
differ only slightly, usually by a few hundredths of a percent for total energies, a few 
tenths of a percent for ehomo and a few percent for the multipole moments. The 
largest difference between the enoMO values occurs for Be2, where they differ by 
3%. For N2, the energetic order of the IvTu and Sug orbital is reversed in all DFT 
approaches as compared to the HF result, which corresponds to the experimentally 
observed order of the outer valence ionization potentials [26|. As far as the multipole 
moments are concerned, the largest discrepancy between the x-only KLI and HF 
approaches occurs for the total hexadecapole moment of N2, where the results differ 
by 8.8%. In this case, the Slater approximation gives a value differing only by 3.7% 
from the HF one. The 1/r expectation values obtained with the HF and x-only KLI 
method are almost identical, differing by only a few hundreths of a precent with 
the exception of the one for the Hydrogen nucleus in FH, where the difference is an 
order of magnitude larger. In this case, both the Slater as well as the x-only LDA 
approximations give values closer to the HF results. 

The Slater method gives - with a few exceptions mentioned above - values for 
Etot. ^HOMO.the multipole moments and 1/r expectation values which differ to 
a larger extent from both the KLI and HF results than the latter from each other. 
From the energy eigenvalues of the inner orbitals, on the other hand, it is obvious 
that the Slater exchange potential V-^^{r) is deeper than the one obtained in the 
KLI method, giving results closer the HF ones. 

Finally, the x-only LDA results differ more strongly from the other methods, 
yielding much higher total energies. The difference is most pronounced for the 
values of shomo. which are roughly twice as large as the ones from any of the 
other methods. This is due to the wrong exponential decay of V^J^^{r) for large 
r. 

We point out that the bulk part of the difference between the x-only KLI and 
the HF results are not due to the KLI approximation, but to the different nature 
of the HF and the DFT approaches. This is an established fact for atomic systems 
p, 10, 11| and we see no reason why it should not hold for molecular systems as 



well. We mention that the difference between the HF and the exact x-only DFT 
results also implies that the exact quantum chemical correlation energy and the 
exact DFT correlation energy are not identical [^] . 
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4 Correlation contributions to the OEP 



The inclusion of correlation effects into the OEP scheme is straightforward, as 
anticipated by the indices xc in section ^, once an explicit functional for Ec [{^jr}] 
has been specified. It has been shown |15, ^] that the orbital-dependent Colle- 
Salvetti functional 29| is well suited for this purpose. It yields excellent results 
for atoms, surpassing the accurracy of conventional Kohn-Sham calculations. In 



this approximation, Ec is given by ||15| 



Ec{W,r}\ = / 7(r)C(r) 



Y,p„{v)Y,\V^U^)\' - -|Vp(r) 



-a j cfr 7(r 



7]{r) 



where 



7(r) 
r/(r) 

e(r) 



1 

1 + dp{r) 3 



Pir) 



and 



a = 0.04918, 
c = 0.2533, 



r]{r) 

b = 0.132, 
d = 0.349. 



(32) 

(33) 
(34) 

(35) 



4.1 Two-Electron Systems 



In order to study the correlation contributions more thoroughly, we first concentrate 
on two-electron atoms for two reasons. First of all, as pointed out above, the 
solution of the full OEP integral equation for these systems is identical to the one 
obtained from the KLI-scheme. As the exchange energy functional is also known 
exactly, c.f. equation (20), the only error made is due to the approximation for Ec- 
Secondly there exist practically exact solutions |^ of the two-particle Schrodinger 
equation so that various DFT-related quantities of interest can be compared with 
exact results. 

In Table IC we show the total absolute ground-state energies of the atoms 
isoelectronic with helium. The first column, denoted by KLICS, displays the results 
obtained with the above described method, including the Colle-Salvetti functional 
for Ec into the OEP scheme. The next two columns show results obtained with 
the conventional Kohn-Sham method for comparison. BLYP denotes the use of 
the exchange-energy functional by Becke |32] combined with the correlation energy 
functional by Lee, Yang and Parr 133], whereas the third column headed PW91 



refers to the generalized gradient approximation by Perdew and Wang [g4|. The 
exact nonrelativistic results in the last column are taken from [30|. Note that there 
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KLICS 


BLYP 


PW91 


exact 




0.5189 






0.5278 


He 


2.9033 


2.9071 


2.9000 


2.9037 


Li+ 


7.2803 


7.2794 


7.2676 


7.2799 


Be2+ 


13.6556 


13.6500 


13.6340 


13.6556 


B3+ 


22.0301 


22.0200 


21.9996 


22.0310 




32.4045 


32.3896 


32.3649 


32.4062 




44.7788 


44.7592 


44.7299 


44.7814 


06+ 


59.1531 


59.1286 


59.0948 


59.1566 




75.5274 


75.4981 


75.4595 


75.5317 


Ne8+ 


93.9017 


93.8675 


93.8241 


93.9068 


Na9+ 


114.2761 


114.2369 


114.1886 


114.2819 


Mgio+ 


136.6505 


136.6064 


136.5531 


136.6569 




161.0250 


160.9758 


160.9175 


161.0320 


Sii2+ 


187.3995 


187.3453 


187.2819 


187.4070 


pl3+ 


215.7740 


215.7147 


215.6462 


215.7821 


S14+ 


246.1485 


246.0842 


246.0105 


246.1571 


C|15+ 


278.5231 


278.4536 


278.3748 


278.5322 


Ari6+ 


312.8977 


312.8231 


312.7390 


312.9072 


K17+ 


349.2723 


349.1926 


349.1032 


349.2822 


Cai8+ 


387.6470 


387.5620 


387.4674 


387.6572 


A 


0.0053 


0.0450 


0.0943 





Table 10: Total absolute ground-state energies for the Helium isoelectronic 
series from various self-consistent calculations. A denotes the mean absolute 
deviation from the exact values from [^] . All numbers in atomic units. 

is no convergence for neagtive ions in the conventional Kohn Sham method. All of 
our calculations have been performed with a basis-set-free, fully numerical atomic 
code which solves the radial Schrodinger equation (||) by the Numerov-method as 
described in ||3^. The angular parts are treated analytically. 

In Figure ]l[ we have plotted the errors -E^t^^ ~ ^tot'^^ °f the numbers shown 



in Table IC. It is obvious that the KLICS scheme gives superior results. The mean 
absolute error, denoted by A, is smaller by an order of magnitude for the KLICS 
results as compared to the two conventional Kohn Sham approaches. This comes 
as no surprise, as the exchange part is treated exactly in the OEP method, whereas 
only approximative functionals can be used in the Kohn-Shame scheme. As may 
be read off Table |l^, where we show the exchange and correlation contributions to 
the total energy seperately for systems where exact values are available, an error 
cancellation occurs in the BLYP and PW91 approaches - the exchange energies 
being too large and the correlation energies being to small in magnitude. The 
KLICS results for these two quantities are clearly much better. However, for the 
highly charged two-electron ions the quality of the results decreases substantially 
in all approaches. For these systems the LYP-functional appears to perform best. 
In order to assess the quality of the xc potentials resulting from various approx- 



12 



o 

CO 



LU 
I 

H 
LL 
C 

LU 



0.20 
0.16 
0.12 
0.08 
0.04 
0.00 
-0.04 



+ KLI-CS 
□ BLYP 
A PW91 



^ S + + + + + + + ^ ^ 



+ + + + 



_I I I I I I L_ 



_I I L_ 



Figure 1: Energy differences corresponding to Table 10 



2 4 6 8 10 12 14 16 18 20 

Atomic Number 





element 


KLICS 


BLYP 


PW91 


exact 




H- 


0.4053 






0.3809 




He 


1.0275 


1.0183 


1.0095 


1.0246 


—Ex 


Be2+ 


2.2768 


2.2573 


2.2367 


2.2766 




Ne8+ 


6.0272 


5.9749 


5.9189 


6.0275 






49.7779 


49.3412 


48.8806 


49.7779 




H- 


0.0312 






0.0420 




He 


0.0416 


0.0437 


0.0450 


0.0421 




Be2+ 


0.0442 


0.0493 


0.0530 


0.0443 




Ne8+ 


0.0406 


0.0504 


0.0615 


0.0457 




Hg78+ 


0.0276 


0.0506 


0.0805 


0.0465 



Table 11: Exchange and correlation energies from various approximations. 



Exact values from |31|. All values in atomic units. 
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imate functionals, it is informative to look at the highest occupied orbital energy of 
the system. In an exact implementation of DFT this value should be equal to the 
ionization potential of the system. Therefore, the resulting values from approxi- 
mate schemes are an indication of the quality of the corresponding xc potential. 
From Table [l^, where we have listed these numbers for various self-consistent ap- 
proximations together with the exact ones, it is obvious that the KLICS scheme 
performs much better than the conventional Kohn-Sham schemes. The difference 
is less pronounced for the highly charged ions as the nuclear potential becomes 
the dominant term in the Kohn-Sham potential (|^). A glance at the second col- 
umn, in which we give the corresponding values from an x-only KLI calculation, 
shows, however, that the reason for the superior quality is due to the inclusion of 
the exact exchange in the KLI scheme, which results in the correct — 1/r asymp- 
totic behaviour of the xc potential. In fact, adding the Colle-Salvetti formula for 
the correlation energy slightly worsens the results, as may be seen by comparing 
the second and third columns: The correlation contribution lowers the already too 
small values from the x-only calculations for the highest occupied orbital energy 
even more. 





KLI 


KLICS 


BLYP 


PW91 


exact 




x-only 


xc 


xc 


xc 




He 


0.9180 


0.9446 


0.5849 


0.5833 


0.9037 


Li+ 


2.7924 


2.8227 


2.2312 


2.2269 


2.7799 


Be2+ 


5.6671 


5.6992 


4.8760 


4.8701 


5.6556 


B3+ 


9.5420 


9.5751 


8.5201 


8.5129 


9.5310 


(-4+ 


14.4169 


14.4507 


13.1638 


13.1554 


14.4062 




20.2918 


20.3261 


18.8072 


18.7978 


20.2814 


06+ 


27.1668 


27.2014 


25.4504 


25.4401 


27.1566 




35.0418 


35.0766 


33.0935 


33.0823 


35.0317 


Ne8+ 


43.9167 


43.9517 


41.7366 


41.7245 


43.9068 


Na9+ 


53.7917 


53.8269 


51.3796 


51.3666 


53.7819 


Mgio+ 


64.6667 


64.7020 


62.0225 


62.0086 


64.6569 


A|ii+ 


76.5417 


76.5770 


73.6654 


73.6506 


76.5320 


Sil2+ 


89.4167 


89.4521 


86.3083 


86.2926 


89.4071 


pl3-H 


103.2917 


103.3272 


99.9511 


99.9345 


103.2821 


514+ 


118.1666 


118.2022 


114.5939 


114.5764 


118.1571 


C|i5+ 


134.0416 


134.0773 


130.2367 


130.2183 


134.0322 




150.9166 


150.9523 


146.8795 


146.8602 


150.9072 


K17+ 


168.7916 


168.8273 


164.5223 


164.5021 


168.7822 




187.6666 


187.7024 


183.1650 


183.1439 


187.6572 



Table 12: Absolute highest occupied orbital energies from various self- 
consistent calculations. Exact values calculated from ||3^. All values in atomic 
units. 

The error in the correlation potential responsible for this behaviour is clearly 
visible from Figure |2| where we plot the exact [31| and various self-consistent corre- 
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Figure 2: Comparison of the exact and self-consistently calculated correlation 
potentials of helium 

lation potentials. The potential obtained with the Colle-Salvetti functional within 
the KLICS scheme shows the same deficiencies as the conventional density function- 
als: Instead of a maximum with positive values of the potential, the approximations 
possess one or even two minima and spurious divergences occur at the origin, which 
may be traced back to gradients of the density and of the one-particle orbitals occur- 
ing in the various correlation energy functional. The need for further improvement 
of the correlation energy functional in this respect is obvious. 



4.2 Beryllium and Neon Isoelectronic Series 

For further analysis we have calculated the total ground state energies of positive 



ions isoelectronic with Beryllium (shown in Table 13) and Neon (shown in Table 
[l4| ). Again, we compare the various DFT methods with exact data from Ref. 

The errors are plotted in Figures |3| and |^, respectively. The data for both 
series show the same trends: The overall best results are obtained with the KLICS 
scheme, where the absolute mean deviation A from the exact values is smallest. 
The BLYP scheme is only slightly worse, but the PW91 functional gives errors 
which are roughly twice as large as compared to the other DFT approaches. From 
the plots in Figures |3| and ^ it is obvious that these statements hold for most ions 
individually. 

There are two other trends to be noted in these results. First we point out 
that although the absolute errors rise within the isoelectronic series as the atomic 
number increases, the percentage errors remain almost constant. And secondly, the 
mean absolute error is smaller by almost an order of magnitude for the ten-electron 
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KLICS 


BLYP 


PW91 


exact 


Be 


14.6651 


14.6615 


14.6479 


14.6674 


B+ 


24.3427 


24.3366 


24.3160 


24.3489 


0+ 


36.5224 


36.5143 


36.4881 


36.5349 




51.2025 


51.1927 


51.1618 


51.2228 


04+ 


68.3825 


68.3713 


68.3362 


68.4117 


F5+ 


88.0624 


88.0499 


88.0110 


88.1011 


Ne6+ 


110.2420 


110.2285 


110.1859 


110.2909 




134.9216 


134.9071 


134.8610 


134.9809 




162.1010 


162.0857 


162.0361 


162.1710 


A|9+ 


191.7803 


191.7642 


191.7113 


191.8613 


Siio+ 


223.9595 


223.9427 


223.8864 


224.0516 


pii+ 


258.6387 


258.6212 


258.5616 


258.7420 


512+ 


295.8178 


295.7996 


295.7367 


295.9324 


C|13+ 


335.4968 


335.4781 


335.4119 


335.6229 


Ari4+ 


377.6758 


377.6566 


377.5870 


377.8134 


K15+ 


422.3548 


422.3350 


422.2621 


422.5040 




469.5338 


469.5134 


469.4372 


469.6946 




519.2127 


519.1919 


519.1122 


519.3851 


Tii8+ 


571.3917 


571.3703 


571.2873 


571.5757 


V19+ 


626.0706 


626.0487 


625.9623 


626.2663 




683.2495 


683.2271 


683.1373 


683.4570 




742.9284 


742.9056 


742.8123 


743.1476 


Fe22+ 


805.1072 


805.0840 


804.9873 


805.3382 


Co23+ 


869.7861 


869.7624 


869.6623 


870.0289 




936.9650 


936.9408 


936.8373 


937.2195 


A 


0.1180 


0.1352 


0.1973 





Table 13: Total absolute ground-state energies for the Beryllium isoelectronic 
series from various self-consistent calculations. A denotes the mean absolute 
deviation from the exact values from [^] . All numbers in atomic units. 

series compared to the four-electron series. 

The ionization potentials from the various approaches as calculated from the 



highest occupied Kohn-Sham orbitals are shown in Tables 15 and 16 for the four- 
and ten-electron series, respectively. The exact nonrelativistic values have been 
calculated from the data given in |Q. Due to the correct asymptotic behaviour of 
the xc-potential for large r within the OEP scheme it comes as no surprise that the 
KLICS data are superior to the conventional Kohn-Sham approach. The effect of 
the correlation potential within the OEP scheme is - like in the two-electron case - 
a lowering of the energy eigenvalue of the highest occupied orbital, as may be seen 
from a comparison of the second and third columns showing the OEP results in 
x-only approximation and with inclusion of correlation in the form of Colle-Salvetti 
in the KLI-scheme, respectively. As opposed to the Helium and Neon isoelectronic 
series, this effect improves the quality of the results in the Beryllium isoelectronic 
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KLICS 


BLYP 


PW91 


exact 


Ne 


128.9202 


128.9730 


128.9466 


128.9376 


Na+ 


162.0645 


162.0956 


162.0668 


162.0659 


Mg2+ 


199.2291 


199.2448 


199.2136 


199.2204 


A|3+ 


240.4071 


240.4102 


240.3768 


240.3914 


SU 


285.5945 


285.5867 


285.5509 


285.5738 


p5+ 


334.7888 


334.7712 


334.7331 


334.7642 


56+ 


387.9885 


387.9616 


387.9212 


387.9608 


cr+ 


445.1922 


445.1567 


445.1138 


445.1622 




506.3993 


506.3554 


506.3101 


506.3673 


K9+ 


571.6091 


571.5570 


571.5092 


571.5754 




640.8211 


640.7610 


640.7107 


640.7861 


Sc"+ 


714.0350 


713.9671 


713.9141 


713.9988 


Tii2+ 


791.2504 


791.1748 


791.1191 


791.2132 


Vl3+ 


872.4671 


872.3839 


872.3255 


872.4291 


Cri4+ 


957.6850 


957.5942 


957.5331 


957.6463 




1046.9039 


1046.8056 


1046.7417 


1046.8646 




1140.1237 


1140.0179 


1139.9511 


1140.0838 




1237.3442 


1237.2309 


1237.1613 


1237.3039 




1338.5654 


1338.4447 


1338.3722 


1338.5247 


A 


0.0293 


0.0334 


0.0694 





Table 14: Total absolute ground-state energies for the Neon isoelectronic series 
from various self-consistent calculations. A denotes the mean absolute deviation 



from the exact values from |36| . All numbers in atomic units. 

series. We mention that the ionization potentials are in much better agreement 

with the exact results if they are calculated as ground-state energy differences [p^]. 
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Figure 3: Energy differences corresponding to Table 13. 
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Figure 4: Energy differences corresponding to Table 14 
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KLI 


KLICS 


BLYP 


PW91 


exact 




x-only 




xc 


xc 


xc 




Be 


0.3089 





3294 


0.2009 


0.2072 


0.3426 


B+ 


0.8732 





8992 


0.7129 


0.7185 


0.9243 




1.6933 


1 


7226 


1.4804 


1.4856 


1.7594 




2.7659 


2 


7975 


2.5000 


2.5049 


2.8459 


^ A 1 

04+ 


4.0898 


4 


1231 


3.7706 


3.7754 


4.1832 


p5+ 


5.6644 


5 


6991 


5.2918 


5.2964 


5.7708 


Ne^+ 


7.4896 


7 


5253 


7.0633 


7.0678 


7.6087 


K 1 7-1- 

Na^+ 


9.5652 


9 


6017 


9.0850 


9.0894 


9.6967 




11.8910 


11 


9281 


11.3569 


11.3611 


12.0348 


A |9+ 


14.4669 


14 


5047 


1 0~700 

13.8788 


13.8829 


14.6230 


Siio+ 


17.2930 


17 


3312 


16.6508 


16.6548 


17.4613 


pii+ 


20.3692 


20 


4078 


19.6729 


19.6768 


20.5496 




23.6955 


23 


7345 


22.9450 


22.9487 


23.8880 


C|13+ 


27.2718 


27 


3111 


26.4671 


26.4707 


27.4763 




31.0982 


31 


1378 


30.2393 


30.2427 


31.3147 


K15+ 


35.1747 


35 


2145 


34.2615 


34.2647 


35.4031 




39.5012 


39 


5412 


38.5336 


38.5367 


39.7416 




44.0777 


44 


1179 


43.0558 


43.0587 


44.3300 


Tii8+ 


48.9043 


48 


9447 


47.8280 


47.8307 


49.1684 


Vl9+ 


53.9808 


54 


0214 


52.8502 


52.8527 


54.2569 


Cr20+ 


59.3074 


59 


3481 


58.1224 


58.1247 


59.5954 


Mn2i+ 


64.8840 


64 


9249 


63.6447 


63.6467 


65.1838 


Fe22+ 


70.7107 


70 


7516 


69.4169 


69.4187 


71.0223 


Co23 + 


76.7873 


76 


8284 


75.4391 


75.4408 


77.1108 


Ni24+ 


83.1140 


83 


1551 


81.7113 


81.7128 


83.4493 



Table 15: Ionization potentials from highest occupied Kohn-Sham orbital en- 
ergies for the Beryllium isoelectronic series from various self-consistent calcula- 
tions. Exact nonrelativistic values calculated from |3£] . All numbers in atomic 
units. 
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KLI 


KLICS 


BLYP 


PW91 


exact 




x-only 


xc 


xc 


xc 




Ne 


0.8494 


0.8841 


0.4914 


0.4942 


0.7945 


Na+ 


1.7959 


1.8340 


1.3377 


1.3416 


1.7410 


Mg2+ 


3.0047 


3.0450 


2.4531 


2.4579 


2.9499 


A|3+ 


4.4706 


4.5125 


3.8285 


3.8339 


4.4161 


S\^+ 


6.1912 


6.2343 


5.4601 


5.4661 


6.1371 


p5+ 


8.1651 


8.2091 


7.3458 


7.3524 


8.1112 


£6+ 


10.3914 


10.4362 


9.4846 


9.4917 


10.3378 




12.8696 


12.9150 


11.8757 


11.8832 


12.8162 




15.5992 


15.6451 


14.5185 


14.5264 


15.5460 


K9+ 


18.5800 


18.6263 


17.4126 


17.4209 


18.5270 




21.8118 


21.8585 


20.5579 


20.5665 


21.7589 


Scii+ 


25.2943 


25.3413 


23.9541 


23.9630 


25.2416 


Tii2+ 


29.0274 


29.0747 


27.6010 


27.6102 


28.9748 


V13+ 


33.0112 


33.0587 


31.4986 


31.5080 


32.9586 




37.2453 


37.2931 


35.6467 


35.6563 


37.1929 




41.7299 


41.7779 


40.0453 


40.0551 


41.6776 


Fei6+ 


46.4649 


46.5130 


44.6943 


44.7042 


46.4126 




51.4501 


51.4984 


49.5936 


49.6037 


51.3979 




56.6856 


56.7341 


54.7432 


54.7535 


56.6335 



Table 16: Ionization potentials from highest occupied Kohn-Sham orbital en- 
ergies for the Neon isoelectronic series from various self-consistent calculations. 
Exact nonrelativistic values calculated from |36| . All numbers in atomic units. 
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5 Conclusions 



Our calculations for molecular systems reveal that the KLI approach is also fea- 
sible for more complex systems and gives results of similar quality as for atoms. 
We expect that an inclusion of correlation effects will result in a highly accurate 
DFT scheme. The studies of correlation contributions to atomic systems show 
that further improvement of the presently available correlation-energy functionals 
is necessary. 
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